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[1] Mathematik unterrichten mit DERIVE, Bernhard Kutzler 
Addison - Wesley, 1995, ISBN 3 89319 860 1 

[2] Mathematisches Praktikum mit DERIVE, A. Garcia u.a., iibersetzt von L. Klingen 
Addison - Wesley, 1995, ISBN 3 89319 857 1 

[3] Nuevas Tecnologias y Ensenanza de Matematicas, A. Garcia, A. Martinez, R. Minano 
Editorial Sintesis, Madrid, 1995, ISBN 84 7738 283 2 

[4] Nuevas Tecnologias en Geometria, E. Roanes Macias & E. Roanes Lozano 
Editorial Complutense, Madrid, 1994, ISBN 84 7491 531 7 

(This is not a DERIVE book but it is a very interesting book written by father and son Roanes for 
everybody liking geometry. Muchas gracias for your book, Eugenio.) 

[5] Handleiding DERIVE, P.E. J.M. Gondrie & G.A.T.M. van Alst 
Academic Service, Schoonhoven, NL, 1994, ISBN 90 395 0080 0 

[6] Toegepaste computer algebra met DERIVE, C. de Ritter & A.K.H. Overeem-Loohuis 
Academic Service, Schoonhoven, NL, 1994, ISBN 90 395 0194 7 

[7] Toegepaste wiskunde met computeralgebra, M. Kamminga-van Hulsen, P. Gondrie & 

G. van Alst, (140 toepassingsgerichte opgaven, uitgewerkt in DERIVE and Maple) 

Academic Service, Schoonhoven, NL, 1994, ISBN 90 6233 956 3 

One of the results of the German DUG Meeting in Nuremberg was the idea to establish an exchange 
for DERIVE teaching materials in the DNL. If you have materials to offer then please write, fax or 
call. I will announce it in the DNL, and if there are other members who are interested in the papers or 
diskettes they should organize contact with each others. 

I can offer: The Binomial Theorem, GCD & LCM, The System of Coordinates (in English and in 
German). 


Modellversuch 

Entwicklung und Erprobung von Materialien 
zum Einsatz von Computeralgebra Systemen 
im Mathematikunterricht der Sekundarstufe II 

Das jeweilige Computeralgebra System soil - wie heute der gewohnliche Taschenrechner - als standig prasentes 
Werkzeug beim Mathematik Lernen Verwendung fmden. Die Materialien (Medienpakete aus Buch & Disketten) 
werden die Funktion eines vollstandigen Mathematik-Lehrwerks iibernehmen, d.h. den gesamten Mathematik- 
unterricht neuer Form unterstiitzen. Sie sollen ferner innovativ auf die Lehrplanentwicklung wirken. 

Solange allerdings die herkommlichen Lehrplane noch in Kraft sind, miissen in mannigfacher Weise Kompro- 
misse eingegangen werden. 

Am Modellversuch sind Lehrkrafte ausgewahlter Celler Gymnasien und insbesondere Mitglieder des Staatlichen 
Studienseminars Celle sowie weitere Versuchsschulen beteiligt. Der erste Materialienband zur Analysis 1 zu- 
nachst fur das Computeralgebra System DERIVE wird zu Beginn des Schuljahres 1995/96 vorliegen (Erpro- 
bungsfassung). Er wird spater in einem namhaften Schulbuchverlag publiziert werden. 

Eine - noch unvollstandige - Vorfassung (100 Seiten, Ergebnis der Vorlaufstudien) kann ab sofort von interes- 
sierten Einzelpersonen oder Schulen zum Selbstkostenpreis (10.- DM in Briefmarken) bezogen werden. 

Kontakt: StD Rudeger Baumann, Gymnasium Emestinum, Burgstr. 21, 29221 Celle 
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Liebe DUG-Mitglieder. 

Nach einem arbeitsreichen Friihjahr freuen wir uns 
alle auf den verdienten Urlaub. Wie Sie aus friihe- 
ren Informationen wissen, haben um die Osterzeit 
in Deutschland zwei wichtige Veranstaltungen 
stattgefunden. Die - auch fur DERIVE sehr erfolg- 
reiche - MNU-Hauptversammlung in Numberg gab 
den Rahmen fur das 2. Deutsche DUG- Treffen. 
Herzlichen Dank alien Besuchern des interessanten 
gemeinsamen Nachmittags. Ein Ergebnis ist die 
Einrichtung einer DERIVE-Materialien-Borse im 
DNL. Naheres finden Sie auf der Informationsseite. 
Lieber Wolfgang Propper herzlichen Dank fur die 
Organisation des Treffens. 

Die DERIVE Days Dusseldorf waren von ca. 250 
interessierten Lehrem aus Deutschland. Osterreich, 
Holland und England sehr gut besucht. An dieser 
Stelle herzlichen Dank an Barbel Barzel und Leo 
Klingen fur die professionelle Abwicklung. Ich 
habe mich besonders gefreut, so viele von Ihnen 
wieder oder zum ersten Mai zu treffen. 

Auf Seite 26 finden Sie zwei Produkte aus meinem 
Grafik- Workshop. Einer der nachsten DNLs wird 
hauptsachlich der Geometrie gewidmet sein. 

Das nachste groBe Ereignis wird die DERIVE- 
Konferenz in Hawaii im August sein, zu der sich 37 
Teilnehmer aus 12 Nationen (darunter die VAR 
und Malaysia) angemeldet haben. Ich werde im 
Herbst dariiber berichten. 

Sie finden in diesem Heft einen Beitrag mit etwas 
Hoherer Mathematik iiber die Besselfunktionen. 
Aus Platzmangel muss dieses Mai auf die nachste 
Folge des Kurvenlexikons verzichtet werden, aber 
die Kardiode und ein Beitrag von Peter Baum zum 
Kurvenlexikon liegen schon bereit. 40 Seiten Urn- 
fang sind offensichtlich noch zu wenig. 

Wir legen dieser Nummer eine Quittung fur Ihren 
Mitgliedsbeitrag fur 1995 bei. Bitte verwechseln 
Sie das nicht mit einer Rechnung oder Zahlungs- 
aufforderung - auBer es steht "ERINNERUNG" 
darauf 

AbschlieBend wiinsche ich Ihnen einen schonen 
Sommer und freue mich schon auf die Zusammen- 
stellung des nachsten DNL. 



Dear DUG-Members , 

After a busy spring we are looking forward to 
well-deserved holidays. As you will know from 
earlier information two important events took 
place in Germany at Easter time. The MNU- 
meeting in Nuremberg - very successful for 
DERIVE - was the frame for the 2nd German 
DUG-meeting. Many thanks to all participants 
for joining an interesting afternoon. One result 
is the establishment of a DERIVE-Materials- 
Exchange in the DNL. (See the first page 
please). Dear Wolfgang Propper, many thanks 
for making the meeting possible and success- 
ful. 

About 250 teachers from Germany, Austria, 
Holland and England attended the DERIVE 
Days Dusseldorf. Many thanks from this place 
to Barbel Barzel and Leo Klingen for their pro- 
fessional organization. I had the pleasure to 
meet there so many of you again or the first 
time. 

On page 26 you will find two products from my 
Graphic-Workshop. One of the next DNLs will 
be dedicated mainly to geometric items. 

The next big event will be the Hawaii DERIVE 
Conference in August. 37 ladies and gentle- 
men from 12 nations (UAR and Malaysia 
among them) have registered. I will give a 
report in the next DNL. 

In this issue you will find one contribution of a 
higher mathematical level about Bessel func- 
tions. You will miss the Lexicon of Curves. But 
the Cardiode and Peter Baum's ideas inspired 
by the Lexicon are ready for the next DNL. 
Believe it or not, 40 pages seem to be too few. 

We enclose a receipt for the 1 995 membership 
dues. Please don't misunderstand this receipt 
as an invoice, except you find "REMINDER" on 
it. 

At the end of my letter I wish you a fine sum- 
mer and I am looking forward to producing the 
next DNL for you. 

Sincerely yours 



Ihr Josef Bohm 
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The DER1VE-NE WSLE TIER is the Bulle- 
tin of the DERIVE User Group. It is pub- 
lished at least four times a year with a con- 
tents of 30 pages minimum. The goals of 
the D-N-L are to enable the exchange of 
experiences made with DERIVE as well as 
to create a group to discuss the possibilities 
of new methodical and didactical manners 
in teaching mathematics. 


Editor: Mag. Josef Bohm 
A-3042 Wiirmla 
D'Lust 1 
Austria 

Phone: 43-(0)2275/8207 


Contributions: 

Please send all contributions to the Editor. 
Non-English speakers are encouraged to 
write their contributions in English to rein- 
force the international touch of the D-N-L. 
It must be said, though, that non-English 
articles are very welcome nonetheless. 
Your contributions will be edited but not 
assessed. By submitting articles the author 
gives his consent for reprinting it in the 
DERIVE Newsletter. The more contribu- 
tions you will send, the more lively and 
richer in contents the DERIVE Newsletter 
will be. 


Preview: (Contributions for the next issues): 

Stability of systems of ODEs, Kozubik, SLO 
Prime Iterating Number Generators, Wild, UK 

Graphic Integration, Probability Theory, Linear Programming, Bohm, A 

LOGO in DERIVE, Lechner, A 

DREIECK.MTH, Wadsack, A 

IMP Logo and Misguided Missiles, Sawada, HAWAII 

3D Geometry, Reichel, A 

Parallel- and Central Projection, Bohm, AUS 

Conic Sections and their 3D Visualisation, Fuchs & Bohm, A 

Muller's Method to solve univariate equations a.o., Speck, NZL 

Vector and Vector Indices Sorting, Biryukov, RUS 

Clear Function Parameters Representation, Biryukov, RUS 

Algebra at A-Level, Goldstein, UK 

Tilgung fremderregter Schwingungen, Klingen, GER 

Utility for Complex Dynamic Systems, Lechner, A 

Some Improvements on the Resolution of ODEs, Fuster, E 

Lexicon of Curves (7) - the Cardiod, Weth, GER 

Bezier Curves with DERIVE, Scheu, GER 

and 

Setif, FRA; Vermeylen, Belgium; Lymer, FRA; Leinbach, USA, Aue, GER; 
Wiesenbauer, A; Keunecke, GER; Roeloffs, NL; Baum, GER 


Impressum: 

Medieninhaber: DERIVE User Group, A-3042 Wurmla, D'Lust 1, AUSTRIA 
Richtung: Fachzeitschrift 
Herausgeber: Mag.Josef Bohm 
Herstellung: Selbstverlag 
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G P Speck, Wanganui, New Zealand 

I have received DNL#17 which contains several articles of interest to me. The DERIVE Newsletters 
certainly serve a useful purpose. 

Though elementary, the issue of solving a simultaneous system of linear equations with a number of 
characters too great to Author (edit) all at once seems to bother virtually all students when first en- 
countered. Thus, if this matter has not been raised in previous NEWSLETTERS, I offer the enclosed 
note. 

The following ’’real life” example may be of interest to DUG readers. 

In considering VENN diagrams showing probabilities for events, 32 simultaneous linear equations in 
32 unknown arise in a natural way in considering a Venn diagram for 5 (completely) independent 
events. Such a system of equations is given in the DERIVE printout at the end of this discussion. The 
number of characters is too great to Author (Enter) all at once. However, as usually the case with 
DERIVE there are ways to circumvent a difficulty. In this case one approach is to use APPEND after 
entering the system of equations piecemeal over several lines. 

(I don't give a copy of the printout because it's too big. I try to explain, editor.) 

Name the unknown xOl, x02, ..., x32. This provides an ordering from the 1 st through the 32 nd final 
solution matrix. Then enter in 4 lines 4 times 8 equations, giving (Set Input Mode Word): 

#1: [x24 + xl9 + ...] 

#2: [xll + xl2 + ...] 

#3: [x20 + x25 + ...] 

#4: [xl6 + x22 + ...] 

#5: [#1, #2, #3, #4] 

Then Simplify followed by Solving the system. 

Kurt Schmidt Koln, Germany 

Gestern erhielt ich das Handbuch fur DERIVE 3. Es ist naturlich voller Fehler. Auf Seite 146 ist die 
Formel fur r falsch. 

Als es noch keinen Computer gab und man alles mit dem Rechenschieber berechnen musste, konnte 
man naturlich dieses Problem nicht so elegant behandeln. Wir machten es damals hochst ungeschickt 
und setzteny = ax. Man erhalt dann eine Parameterdarstellung: 

K. Schmidt points out that the manual of DERIVE 3 contains some mistakes. One of them can be 
found on page 1 51 . The formula for r is wrong. In earlier times - without computer support - they did it 
very clumsy and set y = a x, which led to a parameter representation: 

t 4 -i t 4 - r 

-,U- (-10<tf<10). 

1 -t 5 1 -t 5 _ 

The given problem is: Implicit plot of x 5 + x 4 = y 5 + y 4 . Substituting r cos 9 for x and r sin 6 for y, then 
solving for r gives two explicit solutions: 

l-2cos 2 0 -cos 20 

r = or r = . 

cos 5 Q — sin 5 6 cos 5 0 — sin 5 6 


t 4 -\ 
1 -t 5 


See the DERIVE realisation (version 6.1) in the following: 
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5 4 5 4 

#1: x + x = y + y 

5 5 4 45 54 4 

#2: r -COS(e) + r -COS(e) = r -SIN(e) -+- r -SIN(0) 

5 54 45 54 4 

#3: SOLVE (r -C0S(9) + r -COS(e) = r -SIN (9) + r -SIN (9) , r) 

2 

1 - 2-COSC9) 

#4 : r = v r = 0 

5 5 

C0S(9) - SIN(0) 


#5: 

# 6 : 

#7: 

#S : 

#9: 


Trigonometry := Collect 
Tri'gpower := Cosines 
- C0S(2-9) 

r = v r = 0 

5 5 

C0S(9) - 5IN(0) 

5 4 5 4 

SOLUTIONS(x +x = y +y a y = t-x, [x , y]) 

0 0 

3 2 3 2 

t+t+t + 1 t- (t + t + t + 1) 


4 3 2 4 3 2 

t+t+t+t + 1 t+t+t+t + 1. 



implicit plot polar plot parametric plot 

Comparing the quality of the plot shows that 

(1 ) only the implicit plot presents the solution y = x and 

(2) the best quality plot of the curve gives the polar form. 

Additional comment from Kurt Schmidt: 

The formula on page 21 in revised DNL#1 1 (Digital Filters Design) can be simplified. Instead of: 
TRIW(k,n) := (ABS(k+n+l)/2+ABS(k-n-l)/2-ABS(k)) WR(k, n)/(n+l) one can write: 
TRIW(k, n) := (n+l-ABS(k)) WR(k, n)/(n+l). 
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Pam Bishop sent some messages from the email DERIVE News Group which is established at the 
Birmingham University. Many thanks, Pam, and the best to you. 


Here is a question from Alfonso Poblacion together with A. van der Meer's answer: 

If we try to calculate 

INT (COS (SQRT (x) ) ,x) 

DERIVE gives us 

2 *COS (SQRT (x) ) +2*SQRT (x) *SIN (SQRT (x) ) 

that is right. 

Then we want to find the area of the region between the graph of y = cos(sqrt(x)) and the 
x-axis in [0, 4k 2 ]. If we try with 

INT (COS (SQRT (x) ) , x, 0, 4*pi A 2) 

0 

the answer is 0, that is also correct, But if we try 

INT (ABS (COS (SQRT (x) ) ) , x, 0, 4*pi A 2) 

0 

the answer is still 0, and this is absurd. If you try with many other functions which intersect 
the x-axis, like sin(x), cos(x), etc. the answers are perfect. What happens with this function? 
Maybe there will be others. 

Later versions of DERIVE give the correct result: 

4- jt a 2 

r | COSC/x) | dx = 8 -tt 

’ 0 



A. van der Meer's comment was very interesting in 1995. I reactivated DERIVE (version 
3. 14) and reproduced his answer: 

One has to be very cautious in integrating values. DERIVE is clever, but it is only a machine. Let us 
see what happens: 

4-7r A 2 

#1: / |C0SG/x)| dx = 0 

0 

#2: gives the wrong answer. Without the absolute value, the 

primitive is correct. 

#3: / COS(Vx) dx 

#4: 2-COSG/x) + 2-Vx-SIN(7x) 
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#27: H(t) := INT_SUBST 


2 - SIN(t) 


, t, TAN 


2 J 


#28: 


2-73-ATAN 


73 • (COS(t) - 2 • SIN(t) + 1) 
3 • (COS(t) + 1) 


#29: 


2-73 -ATAN 


H(t) := - 


73-(COS(t) - 2 • SIN(t) + 1) 
3 • (COS(t) + 1) 


#30: H(tt) 

#31: ? 

#32: H(2-tt) - H(0) 

#33: 0 


H(t) has a discontinuity at t = n. And of course #32 gives the wrong answer 0. 

How do the Tl-calculators perform? 


|'F 17 »ir|T F£T Y FjT Y FHt ’I 

I v |fl 1 gebra |Ca 1 c lother 

Pr gm I Life 1 ean Up] 

n 

A ™ 1 i. 1 1 J 1 1 1 ■— 1 1 ai IU 

■jg * cosC|x)dx 

■jg 71 |cos(Jx)|dx 

■ nInt(|cos(Jx)| , x , G , 4 ■ ji 

■ S ' JI 

G 

jg 71 |cos(Jx)| dx 

2 ) 25.1327412237 

25. 1327412287 

i 

MAIN RAD EXACT 

FUNC 30/30 




|'F17 | ir|T F 1 

|fllgebra| 

f F3t 'i 

Calc 

f FHT Y FE Y FbT Y h 

Other |Prgn 1 0|C lean Up| 

■ 

■ 

■ J 

cosUxJdx 

JI 2 

, 4 
■■ ^ 

4 71 _ cos(-Jx)dx 

0 JI 2 

4 

i-2 + 4ji + 3ji + 

4 ■ ji 
3 ■ ji + 2 

2 8ji 

ii— 2 + 4*n+ 3 #n+ 21 

MAIN RAD EXACT FUNC 30/30 


|'F 17 »ir|T F£T Y F3 t Y FHt Y FE Y Fb^ T "i 

|fl 1 gebra |CalG|other|Prgro I o|c lean Up| 1 

■ 

JI 2 

^ cos(Jx)dx 
' 0TI 2 


TI - 2 

■ 

^ cos(Jx)dx 
TI 2 

4 ■ ji 

... <cos<JXx>> 

- x - n A 2/4 - 9#n A 2/4> >1 

MAIN RAD EXACT FUNC 30/30 


The next messages are from the Bulletin Board Service. As I've heard in many talks when I 
was round in Europe many of our members appreciate the ideas, the hints and the conversa- 
tions of the BBS. 

Message 3392: From KEITH WILLIAMS to HADUD about DEG TO HMS AND HMS TO 
DEG 

Dear Hadud, I came up with functions which might be good to add to your DERIVE functions. You 
may have already made these functions. They are functions which will convert between degrees to 
hours, minutes, seconds and hms to degrees. These two functions use two other functions which ex- 
tract the integer part of a number and the fractional part of a number. 
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ip (m) : =IF (m<0 , -FLOOR ( ABS (m) ) , FLOOR (m) ) 


fp (m) : =IF (m<0 , -MOD (ABS (m) ) , MOD (m) ) 

deg_2_hms (y) :=ip (fp (y) * 60 ) /100+ip (y) +fp (fp (y) *60) *60/10000 
hms_2_deg(y) : =ip (y) +ip ( fp (y) *100) / 60+100*fp ( 100*y) /3600 


Maybe these functions will better improve DERIVE since DERIVE does not have these functions. 
Keith Williams 


Message 3398: From MONA to PUBLIC about INVERSE FUNCTIONS 

When I graph the inverse function of/{x) = x 3 - 1, 1 only get the portion where x > 0. 

Also, if I try to evaluate a function such as the inverse of J{x) which involves a fractional exponent 
such as 2/3 which should yield real values for negative values of x, I get complex answers. 

Note: On my TI-81 calculator I solved this problem by rewriting x A (2/3) as (x A (l/3)) A 2 but DERIVE 
would not Simplify this. Any help would be appreciated. 


Message 3399: From JERRY GLYNN to MONA about #3398 INVERSE FUNCTIONS 

Derive is complex number oriented. Try (— 8) A (l/3) and Simplify and youTl get a complex number and 
not -2. -2 is -2 + 0-/ and would be plotted on the complex plane off to the left at 180 degrees. De- 
rive's answer is 1 + V3 -i which plots in the first quadrant so qualifies as the principal value. 

Your question is the same as mine when I first used Derive (or Mathematica) and Derive has imple- 
mented at my suggestion a bit of flexibility to satisfy us. Do Manage Real (Options > Mode Settings > 
Branch > Real in Derive 6.10) and then simplify (— 8) A (l/3) again and you will get -2 and your plots 
will give what you want. Keep asking questions. 


1/3 


#1: 

C-0) 

= 1 + j3-i 

#Z : 

Branch : 

- Real 


1/3 


#3: 

(-8) 

= -2 


Derive 6.10 



i Real or Complex Format: Rectangular 


- 8 ) 


1 + ' 


/s-i 


Tl-Nspire (same as TI-92 and Tl-Voyage 200) 


Message 3400: From BOOM-BOOM to MONA about #3398 INVERSE FUNCTIONS 

Try (x A 2) A (l/3) for a TI-81 ’’look alike” graph of x A (2/3). The real problem (no pun intended) is the 
complex domain of DERIVE vs the real domain of TI-81. 

Another interesting result is to graph y = abs(ln x). DERIVE and the TI-85 will plot something in the 
second quadrant; TI-81 does not. What's right? (The absolute value - better known as modulus - of a 
complex number is a real number.) Adios, Larry Gilligan 
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Message 3402: From DPHILLIPS to PUBLIC about LARGEST KNOWN PRIME NUMBER 

I was able to compute 2 A 859433 - 1 with Derive XM. Derive could not compute the number directly, 
perhaps because I have only 4 Megs of RAM so I computed it indirectly. I first calculated 2 A 200000, 
multiplied the result by itself, and then multiplied that result by itself to obtain 2 A 800000. I then mul- 
tiplied that result by 2 A 59433 and subtracted 1. The number took 6 minutes and 21.1 seconds to com- 
pute and half an hour and 137 pages to print out. It had the required number of digits, 258716 and the 
number ended in 1. I just hope that Derive was able to compute all of the digits correctly! I gave the 
printout to my daughter, who is in the 1 1 th grade, and is required to keep a math journal. I hope her 
teacher and class appreciate that the 137 page print out is the largest known prime number. 

Don Phillips 

Message 3403: From VICTOR SANTOS to PUBLIC about INVERSE FUNCTIONS 

After reading Jerry Glynn's reply to Mona about the way Derive handles complex branches, I remem- 
bered I found a Derive bug concerning this same subject. I am using Derive 2.55 so I think this bug 
isn't corrected in Derive yet. 

Everyone knows that (-2) A V3 is ALWAYS COMPLEX, no matter what branch you choose. With 
Derive when you Approximate this expression you get the correct answer if Branch is Principal (you 
get a complex number), but when you choose Branch Real you get an erroneous reply (a real number). 
It's interesting to note that you get correct answers in both cases if you SIMPLIFY the expressions 
instead of using APPROX. Obviously there is something wrong here because the answer should al- 
ways be complex. 

Branch : =Real 

(-2) A SQRT (3) 

; Approx (#1) 

2. 21288-2. 47766*#i 

Branch : =Real 

; Approx (#1) 

3.32199 

Should someone from the Derive team see this, I would like to know the reason for this inconsistency. 
Victor Santos, from Portugal, in Europe. 

Message 3404: From HARALD LANG to VICTOR SANTOS about #3403 

No, I don't agree that this is a bug in DERIVE. On the contrary, it is an example on how carefully 
DERIVE is designed. The reason is the following: (-2/, where q is an irrational number, has infinitely 
many values, and they are distributed as a dense subset of the circle of radius 2 q . I.e., every point on 
that circle has a value of (~2) q arbitrary close to it; in other words, any point on that circle is an arbi- 
trary approximation to (~2) q . When DERIVE approximates, using the real branch, she chooses 2 q , 
which is real, and an approximation to any degree of accuracy. IMHO, that is the most relevant choice. 
— Harald 

Message 3406: From VICTOR SANTOS to HARALD LANG about #3405 

Thanks for the explanation about the way Derive handles infinitely many values of (-2) q , q being irra- 
tional, when we approximate using the Real Branch. It now seems reasonable to me to get a real value 
for APPROX of (-2) ' 3 . 
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Nevertheless, something still is not quite clear to me and I again ask for your comment on the follow- 
ing: (2 + i )^ 3 is ALWAYS COMPLEX, too and its infinitely many values all lie on the circle of radius 
5 A (V3/2). Following your explanation, there is a real approximation to any degree of accuracy to 
(2 + i )^ 3 and that is 5 A (V3/2). But in this case Derive will not APPROXimate to this real value, even 
when we choose the Real Branch. 

I would appreciate your comments on this. Thanks again, Victor Santos. 

Message 3408: From HARALD LANG to VICTOR SANTOS about #3406 

Hmm, no I can't see any good reason for this. I suggest you explicitly address your question to Soft 
Warehouse, so they get alerted. Or does Jerry Glynn (Hi, Jerry, I suppose you read this.) have any 
suggestion as to why DERIVE behaves like this? — Harald 

J 3 

#6: (2 + l) 

#7: Branch := Real 

approximated 

#8: 2.798998717 + 2 .399664980 - 1 

simplified 
-/3/2 

#9: 5 

#10: Branch :z Principal 

J3/2 L-C/3-7T/4 - ^3-ATANCl/3» 

#11: 5 -e 

#12: 2.798998717 + 2 .399664930 - l 

Derive 6.10 


Compare how the behaviour is 
changing as time goes by. 

Message 3411: From JERRY GLYNN to PUBLIC about FORM OF ANSWERS IN DERIVE 

In Derive the derivative of certain expressions is handled correctly but in a form which I believe is a 
mistake in judgment. What do you think? 

d 25 2 117 

#1: (x +(l+x+x) ) 

dx 

233 232 231 230 229 228 

#2: 234- x + 27261-x + 1601496-x + 63225162-x + 1886202630-x + 45338551947 -x + 

227 226 225 224 

914299617 348- x + 15904855609578 • x + 243557328977586-x + 3334319272823175 • x + 


#1: (-2) 

#2: Branch := Principal 

#3: 2.212884986 - 2 .477661129 - L 

#4: Branch := Real 

#5: 2.212884986 - 2 .477661129 - L 

Derive 3.14 


UpH 


(2 + 1 

0 J3 


J3 


J3 

5 2 

■ cosCtarrV \s2) ■ J3) + 5 

2 

J3 


J3 

5 2 

■Gos(tan*a^2)-j3) + 5 

2 


l-' 1 

I! tan-iCl-^ 
2. 79899871744 + 2. 39966493676 ■ i 


^S^Iot r r er lp4 E .IolcU^ UpH 

— TT — 

5 2 ■ Gos(tan*( 1^2) -J3) + 5 2 ■ sin(tarrK 1>"^ 
J3 J3 

1 5 2 ■ cos(tan-V 1- / 2) J3) + 5 2 ■ sinCtairVl^ 
2. 79399371744 + 2. 39966493076 ■ i 

J3 

i [•' ^ ^ j f3 p t ~ 1 J 3 . cj 21 
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If I did this expression by hand I never expand the polynomial but Derive does. 

My answer is 25x 24 + 1 17(1 + x + x 2 ) 116 (1 + 2x) and I am done. I think that Derive should not expand 
when asked to differentiate. I am sure there are other sides to this question which I am not seeing. 
Comments? Is this part of a larger question which we might discuss? 

Message 3412: From HARALD LANG to JERRY GLYNN about #3411 

Hi Jerry, 

I am glad you bring this up, since I have had similar problems, and I have also posted messages about 
it here, but never got any responses. I have also asked SWHH if they could explain a little how DE- 
RIVE ’’thinks”, so one could detect what she is trying to do in similar situations and figure out a way 
to prevent her (but no answer, alas). 

For example: if you split the expression and differentiate x 25 and (1+x+x 2 ) 117 separately, then she is 
doing ok, which you have no doubt discovered. Another way out, which is the way I currently use in 
similar situations, is to replace ”117” by, say ”a”, where ”a” is declared a positive number. I.e., I ask 
DERIVE to differentiate a more general function. She will now come up with the expected answer, 
with a term (,..) fl_1 . Now, after you have substituted back 117 for a , you must not simplify the whole 
expression, since then she starts expanding again; rather you must highlight only the exponent ’’in- 
i’’ (from a - 1 above) and simplify that to get the answer we want. 

I think it is extremely important to know that such things happen, and that there are ways out. I.e., 
when you try to do some calculations, and DERIVE just works and works and eventually runs out of 
memory, you should not give up. It seems to me that DERIVE is sometimes ”over-ambitious” as to 
simplifications, and then the trick is to prevent her from even trying, by making it impossible. In this 
case she cannot expand since a is not a number, as opposed to (...) 117 where she tries, until she 
runs out of memory (I don't have the XM-version.) — Harald 

See how other CAS behave in this case (btw, DERIVE has not changed its behaviour since 1995!) 


As you can see on the right screen shot the TIs (V200 and NSpire as well) also try to expand before 
finding the derivative - or try to expand the derivative - but in case of "too much" calculating and dis- 
playing the result for the handheld's memory, both give up and display the unexpanded result - exact 
what we want the device to do. 



ifl^T FF? fht y f? t Fe^ t "i 

I t f — In 1 gebra |Ca 1 c |dther |Pr gn 1 0 |c 1 ean Up | [ 


+ x + x 2 ) 117 ) 

117-(2-x + l)-(x 2 + x + l) 116 + 25 x 24 
■‘^(* 25+ ( 1 +X + X 2 ) 11 ) 


25 x 24 + 22 -x 21 + 231 x 20 + 1320 x 19 + 52> 

B 

IMHIH 


<Kx A 25+<l+x+x A 2> A l:Ux> 


Now let's look how MAXIMA is displaying the derivative ( MATHEMATICA behaves pretty the same!): 


(%il) x A 25+ (l+x+x. A 2) A 117; 

f 2 V 17 25 

(%ol) fx + x+lj +x 
(%i2 } diff(%, x); 

( 2 V 16 2. 

(%o2 ) 117(2 x+l\x +x+l) + 25 x 
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Message 3413: From JERRY GLYNN to HARALD LANG about #3412 

Thank you for your usual thoughtful answer. Somehow, I think all of this is related to the following 
situation: if I have function{s omething) I may want to do the something first and then do function of 
that result OR I may do better to do function of something first and then either carry out the something 
or not. I suppose it gets out of hand if a number of these are nested so I’ll ignore that for the moment. 
If I understand SWHH correctly (I often don't) they believe that Simplify should automatically do the 
right thing and then DERIVE will continue to deserve the ease of use reputation which it has earned. 

I believe this view needs modifying since in some cases it is impossible for DERIVE to know what is 
best. I, in fact, for my own reasons might like to carry out a calculation in more than one way. My 
only solution is for DERIVE to give the user some control of the situation. 

Maybe two versions of Simplify : 

Insimplify and Outsimplify for inside first or outside first. One of the great strengths of DERIVE is its 
appearance of simplicity. If my suggestion were accepted this could be in danger so I understand if 
SWHH is slow to make such a change. Also, it may just be a bad idea. 

Thanks for your contributions to this forum . . . I’m always interested in your comments. 

Message 3420: From ANDROMEDA to PUBLIC about AVERAGES 

A silly math problem has bugging me for some time now which may be simple for most of you out 
there. Text defines the AVERAGE speed as Total Distance/Total Time (a slope). At the same time it 
describes AVERAGE value as (al+a2+a3+ ...)/n or (\l(b-a)]f{x) dx). For example: 


A 

B 

C(=A/B) 

2.40 

5.10 

0.47 

1.30 

3.80 

0.34 

9.10 

15.00 

0.60 

5.90 

6.30 

0.94 

5.00 

6.30 

0.79 


Total A/Total B = 0.65, while Average of column C = 0.63. Both averages have a different value, the 
question is which one of these is valid and WHY? Anyone? 

Message 3421: From JERRY GLYNN to ANDROMEDA about #3420 / AVERAGES 

I say 0.65 is correct. This is the sum of the A's divided by the sum of the B's. Suppose A is the miles 
you travel and B is the number of minutes to do this travel. Adding up the distance and dividing by the 
time added up gives you the average speed for this time in miles per minute. Why doesn't the second 
procedure work? You will get the right answer by the second procedure IF each time segment is the 
same. That way you are averaging speeds which occur for equal times. Otherwise your average fa- 
vours those with a small time. A good trick is to move to extreme cases. If you go 500 miles per hour 
for 23 hours and 10 miles per hour for one hour will your average be (500+10)/24? If you travel 12 
hours at 500 mph and 12 hours at 10 mph then your average would be (500+10)/24. 

Any help? Keep asking questions and you will learn. 
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AN ILLUSTRATION OF THE CONTINUED FRACTIONS 
METHOD IN DERIVE: THE EVALUATION OF BESSEL 
FUNCTIONS OF FRACTIONAL AND INTEGER ORDER 

J. Cerdan, P. Fernandez de Cordoba, D. Ginestar and Yu. L. Ratis 


1 Introduction 

In this paper a double objective is presented in the framework of the Derive system and its tutorial 
utilization. In section 2 we insist in how the continued fractions representation of a given function can 
be used as an approximation of this function, that generally converges faster than another approxima- 
tion as the Taylor expansion. In sections 3 and 4, we introduce an algorithm to calculate Bessel func- 
tions (BFs ), of integer and fractional orders, based on the continued fractions method. This proposed 
algorithm is very efficient numerically [1]. In fact, the usual numerical methods to calculate BFs take 
into account normalization relations [2]. In this paper we introduce an algorithm and corresponding 
Derive code to evaluate regular and irregular BFs without any re-calculation thxough nor- malization 
relations. Furthermore, the method maintains the stability of each recurrence relation, i.e., we use 
forward recurrence relations for the BFs of the second kind and backward ones [3] for the BFs of the 
first kind. The algorithm uses forward recurrence relations to generate irregular BFs and takes into 
account the continued fraction method to evaluate high order regular BFs. From these values we can 
generate regular BFs applying backward recurrence relations. 

The method to evaluate the continued fractions representation of a function and its application to 
the calculation of BFs has been programmed to be used with Derive, which is an exceptionally easy- 
to-use computer system with simple structures and a great pedagogical value. The code has been 
designed in an interactive way and this allows the user to follow the process step by step. 

2 Evaluation of continued fractions 

In refs [4] and [5] we can find a set of operative definitions about continued fractions. 

A continued fraction is given by 



( 1 ) 


( 2 ) 


In the previous expressions the a's and b's can themselves be functions of x. We consider, as an ex- 
ample [4], the following continued fraction representation of the tangent function 


tan(x) 


X x 2 x 2 x 2 

1 ^ 3 ^ 5 ^ 7 ^' 


( 3 ) 


i.e. for this case 
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1. Derive program to evaluate the continued fraction representation of the 
tangent function (see Appendix A). 

The code is organized as follows: 

- Define the a's and b ' s in expression (4). We evaluate the values of a n with the function 

A(N_, z) and the values of b n with the B(N_) function. 

- Define the a's and jB's in expression (5) using the recurrence relations (6) and (7). We 

evaluate the values of a x and jS x with the functions ALPHA(L_, z) and BETA(L_, z). 

a 

- We store the values of f k =—, (k = 0, 1,2, ..., TV) in the H(N,z) vector. 


(/ left all files in the form of 1995 unalterred as much as possible. The first file appears in its 
"DERIVE for DOS" form. Josef.) 

"Appendix A" 

A (N_, z) : =IF ( N_< 1 , 0, IF (N_=l, z, -z A 2) ) 

B ( N ) :=IF (N_<=0, 0, 2*N_-1) 

ALPHA (L_, z) : =IF (L_=-l, 1, IF (L_=0 , B (L_) , IF (L_>=1 , B (L_) * ALPHA (L_-l, z) + 
ALPHA (L_, z) * ALPHA (L_-2 , z) ) ) ) 

BETA (L_, z) : =IF (L_=-l, 0, IF (L_=0, 1, IF (L_>=1,B (L_) *BETA (L_-l, z) + 

A (L_, z ) *BETA (L_-2 , z) ) ) ) 

H (N, z) :=VECTOR (ALPHA (k_, z) /BETA ( k_, z) , k_, 0,N) 

See the results in paragraph 6. 

2. Derive program to evaluate the BF(,-) (see Appendix B). 


The code is organized in the following way: 

- Evaluate the set (y„(z), n = 0,1, 2, Nmax} using the forward recurrence relation (8) 

and the known values of yo(z) and y\(z). We evaluate r/r) with the PAS0(L_, z) func- 
tion, and we store the set of calculated values in the YNEW(NMAX.z) vector. In order 
to compare our results, we also evaluate this set using the 
SPHERICAL_BESSEL_Y(k_,z) function defined in the standard BESSEL.MTH utility 
file - wihich is now BesselFunctions.mth - storing the results in the YOLD(NMAX,z) 
vector. 

- Use of the continued fraction method to evaluate H(z ) of expression (9) which we rewrite 

in the form 


Pk 


H(z)=- 


-1 


-1 


-1 


( 20 ) 



' max + -j 


)z ~ 1 + 2(TV> 


' max + f 


'jz + 2 ^JVlTiClX ~h Z + 


and so, 


n > 1 


b () = 0, b n =-(2Nmax + 2n-\); n> 1 


z 


( 21 ) 


D 
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We calculate the values of a n and b n of the expression (21) with the A(N_) and 
B(N_, NMAX , z) functions, and the values of a t and Pi with the ALPHA(L_, NMAX, z) 
and BETA(L_ , NMAX , z) functions. 


The last step is to store the values of fk(k=0, 1,2, ..., N_) in the H(N_, NMAX, z) vec- 
tor. 


By looking at the elements of the H (N_ , NMAX) vector we get the value of H(z): 

H{z) = ELEMENT(H(N_, NMAX, z) , nl im) 

and nl i m is defined as the component of the vector H (N_ , NMAX , z) where the 
changes are small with respect to the previous components. 

- Evaluate /\„ MV i(~) and j Nmax (z). 1T0PM1(NMAX , z) is a function to calculate /y,„ l7V i(-) 

using the equation (10) and evaluate /. v „ MV (z) using the expression (11) with the func- 
tion JT0P(NMAX , z) . 

- Evaluate the set {j„(z), n = 0, 1, 2, ..., Nmax} using the backward recurrence relation 

(12). We evaluate ji(z) with the PASO! (L_, z) function, and we store the set of calcu- 
lated values in the 1 NEW (NMAX , z) vector. In order to compare our results, we also 
evaluate this set using the SPHERICAL_BESSEL_l(k_,z) function, defined in the 
standard BESSEL.MTH utility file, storing the results in the 10LD(NMAX , z) vector. 


#1: LOAD(C:\Programme\TI Education'.Derive 6\Mafh\BesselFunctions .mth) 

Appendix B - BesselFunctions.mth (earlier BESSEL.MTH) are preloaded 
Bessel Functions of fractional order (Spherical Bessel Functions) 

PASO (I . z) 

If L_ = 0 
- CQS(z)/z 
#2: If L_ = 1 

(- SIN(z) - C0S(z)/z)/z 

(2-L_ - l)/z -PAS0(L_ - 1, z) - PAS0(L_ -2, z) 

#3: YNEW(NMAX, z) := VECTOR (PASO (k_, z) , k_, 0, NMAX) 

#4: YQLD(nmax r z) := VECTOR (SPHERICAL_BESSEL_Y(k_ r z) , k_ r 0, nmax) 

Use of the continued fractions method 

#5 : InputMode := Word 

A(N_) := 

If N_ > 0 

#6: -1 

0 

B(N_ ? \mX. z) := 

If N_ = 0 
0 

#7: If N_ > 1 

1/z- [2 ■ (NMAX +■ N_) - 1) 

0 

ALPHA(L_ P NHAXp z) := 

If L_ = -1 
1 

If L_ = 0 

B(L_. NMAX r z) 

If L_ > 1 

B(L_. NMAX, z) ■ALPHA(L_ - l f NHAX ? z) # A(L_) -ALPHA(L_ - 2. NM^X. z) 
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BETA(L_, NHAX, z) 

If I = -1 

0 

#9: If L_ = 0 

1 

If L_ > 1 

B(L_. NHAX, z) ■BETA(L_ - 1, NHAX, z) + A(L_) .BETA(L_ - 2. NHAX, z) 


( ALPHA(k_, NHAX, z) 

#10: H(N_. NHAX, z) VECTOR 

K BETA.(k_, NHAX, z) 

Spherical Bessel functions of the first kind 


, k_, 0, N_J 


JT0PH1CNHAX, z) 

# 11 : 2 


z -((HCnl-im, NHAX, z)) ■ (YNEW(NMAX, z)) - (YNEW(NMAX, z» ) 

nlim NHAX NHAX + 1 


J TOP (NHAX, z) := (H(n1im, NHAX, z)) OTOPH1(NHAX, z) 

#12 : nlim 

pasq] (l_„ nhax. z) := 

If L_ > NHAX 
0 

If L_ = NHAX 
JTOP (NHAX . z) 

#13: If L_ = NP-1AX - 1 

]T0PH1(NHAX, z) 

If L_ > 0 

(2 -L_ +■ 3) /z - PASO] (L_ + 1 P N MAX . z) - PASO] (L_ + 2. NHAX, z) 
0 


#14: ]NEW(NMAX P z) VECTOR (PASO] (j_ P NHAX P z) p j_ P NHAX P 0 P -1) 

#15: ]QLD(NMAX P z) := VECTOR(SPHERICAL_BESSEL_] (k_ P z) P k_ P N HAX , 0 P -1) 


3. Derive program to evaluate the B(s) (see Appendix C). 

The code is organized in the same way as in the SBFs case: 

- Evaluate the set {Y n (z) 9 n = 0, 1,2, ..., Nmax} calculating 7 0 (z) and Y\(z) and using for- 

ward recurrence relation (15). 

To evaluate Y 0 (z) and 7i(z) we use the expressions (13) and (14). These expressions 
are given in terms of infinite series. To get these values we define the vectors 
VECTOR_YO(N_,z) and VECT0R_Y1(N_, z) where we store the different partial sums 
up to order N_. The values 7 0 (z) and 7i(z) are given by the component of these vec- 
tors where the changes are small enough. 

We define the values Y0 = 7 0 (z) and Y1 = 7i(z). As in the previous example, we store 
the set {7„(z), n = 0,1, 2, ..., Nmax } in the YNEW(NMAX, z) vector. 

In the function PAS0(L_, z) we define yi(z). 

- Use of the continued fractions method to evaluate H(z ) of expression (16). 

- Evaluate jumaxi?) and j^max- i( z )* This is done using the expressions (17) and (18). 

- Evaluate the set {/ w (z); n= 1, ..., Nmax). This is done by using the relation (19). 

(Preloading BesselFunctions.mth is not necessary with recent Derive versions. Derive recognizes all 
functions provided in the utility files which are collected in the MATH-folder). 
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#1: LQAD(C:\Programnne\TI Edii cat i on \ Derive 6\Math\BesselFunctions . mth) 


Appendix C 

BesselFunctions.mth (earlier Bessel. mth) preloaded 

Bessel functions of integer order 
Bessel functions of the second kind 

#2 : InputMode := Word 

#3: ]0Cz> :z BESSEL.] (0 r z) 

#4: ]l(z) :z BESSEL.] (1, z) 

#5 : PHI(N_) := -y + HARMDNIC_NUMBER(N_ - 1} 

( 


# 6 : 


# 7 : 


VE CTO R.YQ (NMAX, z) VECTOR 


VE CTO R.Y1 ( N MAX , z) VECTOR 




2 m_ 

LN|| — OO(z) I 

]T k_=0 


PHI(k_ + 1). 


4 J 


m_ r 1, NMAX 


k_! 


2 

JT-Z 




(PHI(k_ + 1) + PHI(k_ - 2)) 


nCz) 


z m_ 

I 

2ir k_=0 


2 ^ 
z 

4 , 


k_ 


k_! ■ (k_ + 1)1 


m_, 1, NMAX 


PASO(L_, z) :z 
If L_ = 0 
yO 

#S: If L = 1 

yl 

(2-L. - 2)/z -PASO(L_ - 1, z) - PASQ(L_ - 2 r z) 

#9: YNEW(NMAX P z) := VECTOR (PASO (k_ ? z) r k_ ? 0 P NMAX) 

Use of the continued fractions method 

A(N_) := 

If N_ > 0 

#10: -1 

0 

B(N_ P N MAX P z) := 

If M_ = 0 
0 

#11: If N_ > 1 

2/z-(NMAX t N_ - 1) 

0 

ALPHA(L_. NMAX, z) := 

If L_ = -1 
1 

#12: If L_ = 0 

B(L_. NMAX, z) 

If L_ > 1 

B(L_ r NMAX P z) ■ALPHA(L_ - l r NMAX r z) + A(L_) -ALPHA(L_ - 2 f NMAX r z) 

BETA(L_p NMAX, z) := 

If L_ = -1 
0 

#13: If L_ = 0 

1 

If L_ > 1 

B(L_„ NMAX, z) -BETA(L_ - l f IMP-lAX P z) + A[L_) -BETA(L_ - 2 P NMAX, z) 


#14: H(N_, NMAX, z) := VECTOR 

Bessel functions of the first kind 

]T0P(NMAX, z) := 


ALPHA(k_, NMAX, z) 
BETA(k_, NMAX, z) 


k_, 0, N_ 


#15: 


ELEMENT (YNEW(NMAX , z) , NMAX) 


ELEMENT (YNEW(NMAX , z) , NMAX + 1) 
ELEMENTCHCnlim, NMAX, z) , nli'm) 
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l 

#16: ]T0PM1(NMAX, z) := ]T0P(NMAX, z) 

ELEMENTCHCnlim, NMAX, z) , nli'm) 

PASO] (L_, NMAX, z) : = 

If L_ > NMAX 
0 

If L_ = NMAX 
] TOP (NMAX, z) 

#17 : If L_ = NMAX - 1 

]T0PM1(NMAX, z) 

If L_ > 0 

(2-L_ + 2)/z- PASO] (L_ + 1, NMAX, z) - PASO] (L_ + 2, NMAX, z) 
0 

#1S: ]NEW(NMAX, z) := VECTOR (PASO] (j_, NMAX, z) , j_, NMAX, 0, -1) 

#15: ]OLD(NMAX, z) := VECTOR (BESSEL_] (j_, z) , j_, NMAX, 0, -1) 


6 Results 

In this section, we have presented some numerical examples of the accuracy of the method. 

1) Results of the program to evaluate tan (z) using a continued fraction representation 

In the previous section we present a algorithm to evaluate the continued fraction representation 
of a given function, and in Appendix A we give as an example, a Derive code to calculate the 
tangent function with the expression (4). In the code, we store the values of^ (see eq. (5)) in the 
H(N , z) vector. The results of tan(z) for the case z = f are the following: 

H(10,pi/4) 

"This the result approximating H (10, pi/4) with DERIVE 6.10" 
[0,0.7853982300,0.9886893656,0.9997873697,1,1,1,1,1,1,1] 

"This was the result in the original MTH-file in 1995 (6 digits)" 

[0,355/452,8479/8576, 4702/4703,1,1,1,1,1,1,1] 
[0,0.785398,0.988689,0.999787,1,1,1,1,1,1,1] 


i.e. we get the right value ( tan(-|) = 1 ) since f 5 . 

2) Results of the evaluation of the BF(s) (See Appendix D). 

Appendix D 

We have Appendix B loaded and we proceed in order to evaluate the 
Spherical Bessel functions of the second kind of all orders below Nmax = 5 
using the proposed algorithm 

#16: YNEW(5, 2) 

approximate #16 

#17: [0.2080734182, -0.3506120042, -0.7339914246, -1.484366557, -4.461291526, -18.59144531] 

In 1 995 an accuracy of 6 digits was used: 

#18: [0.208073, -0.350612, -0.733991, -1.48436, -4.46129, -18.5914] 
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We evaluate the Sph. B.f. of the 2nd kind of all orders below Nmax using 
the utility file BesselFunctions.mth (approximated). 

#19: Y0LD(5, z) 


# 20 : 


COS(z) 


z 


2 2 
COS(z) SIN(z) (z - 3) -COS(z) 3-SIN(z) 3.(2-z - 5)-COS(z) 

2 z 3 2 4 

z z z z 


2 2 4 2 

(z - 15) *SIN(z) 5 • (2 *z - 21) *SIN(z) (z - 45-z + 105).COS(z) 

» — » 
3 4 5 

z z z 


4 2 4 2 

15- (z - 28- z + 63) -COS(z) (z - 105-z + 945)-SIN(z) 


6 5 

z z 

We substitute for z = 2 and approximate 


#21: 


#22: 


COS(2) COS(2) SIN(2) COS(2) 3-SIN(2) 9-COS(2) ll.SIN(2) 59-COS(2) 

2 4 2 ' 8 4 ’ 16 8 ’ 32 

65 *SIN(2) 495 *C0S(2) 541.SIN(2) ' 

16 ’64 32 

[0.2080734182, -0.3506120042, -0.7339914246, -1.484366557, -4.461291526, -18.59144531] 


with 6 digits accuracy as it was in 1995: 

#23: [0.208073, -0.350612, -0.733991, -1.48436, -4.46129, -18.5914] 

We evaluate up to seven terms in the continued fractions expression of 
H(z=2) with Nmax = 5 (first approximated). 

#24: H (7, 5, 2) 

#25: [0, 0.1818181818, 0 . 1870503597 , 0.1871631553, 0.1871649922, 0.1871650156, 0.1871650158, 0.1871650158] 

with 6 digits accuracy as it was in 1995: 

#26: [0, 0.181818, 0.18705, 0.187163, 0.187164, 0.187165, 0.187165, 0.187165] 

We define nlim. 

#27: nlim := 5 

We evaluate the Sph. B.f. of the first kind of order Nmax-1 (approx). 

#28: JT0PM1(5 , 2) 

#29: 0.01407939267 

We evaluate the Sph. B.f. of the first kind of order Nmax (approx). 

#30: IIT0P(5, 2) 

#31: 0.002635169421 

We evaluate the Sph. B.f. of the 1st kind of all orders below Nmax using 
the proposed algorithm (approximated). 

#32: JNEW(5 , 2) 

#33: [0.002635169421, 0 . 01407939267 , 0.06072209763, 0.198447949, 0.4353977749, 0.4546487134] 

Accuracy 6 digits as it was in the original 1995 contribution 

#34: [0.00263516, 0.0140793, 0.0607221, 0.198447, 0.435397, 0.454648] 

We evaluate the Sph. B.f. of the 1st kind of all orders below Nmax using 
the utility file BesselFunctions.mth, substitute for z = 2 and finally approximate. 


#35: J0LD(5, z) 
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#36: 


4 2 4 2 2 

15- (z - 28-z + 63) *SIN(z) (z - 105-z + 945)-C0S(z) 5.(2-z - 21).C0S(z) 


4 2 2 2 2 

(z - 45- z + 105) *SIN(z) (z - 15)-COS(z) 3-(5 - 2-z )-SIN(z) (3 - z )-SIN(z) 

, + , 


3*COS(z) SIN(z) COS(z) SIN(z) 


#37: 


z z 

541«COS(2) 495 *SIN(2) 65-COS(2) 59-SIN(2) ll-COS(2) 9-SIN(2) 3-COS(2) 


32 64 16 

SIN(2) SIN(2) COS(2) SIN(2) 


32 


8 


16 


8 4 2 2 

#38: [0.002635169721, 0 . 01407939275 , 0 . 06072209765 , 0.198447949, 0 .4353977749, 0.4546487134] 

In 6 digits accuracy as it was in 1 995: 

#39: [0.0026324, 0.0140787, 0.0607218, 0.198447, 0.435397, 0.454648] 


In the following table we present some results: 


Function 

Our Algorithm 

BESSEL.MTH(95) 

Reference [5] 

BesselFunctions.mth (08) 

75 ( 2 = 2 ) 

0.002635169421 

0.0026324 

0.0026352 

0.002635169721 

74 ( 2 = 2 ) 

0.01407939267 

0.0140787 

0.014079 

0.01407939275 

73 ( 2 = 2 ) 

0.06072209763 

0.0607218 

0.060722 

0.06072209765 

72 ( 2 = 2 ) 

0.198447949 

0.198447 

0.198447 

0.198447949 

7'i(z=2) 

0.4353977749 

0.435397 

0.435397 

0.4353977749 

7o(z=2) 

0.4546487134 

0.454648 

0.454648 

0.4546487134 


These results are obtained with the same Derive precision (the Derive precision by default). 
(This was 6 digits in 1995 and is 10 digits now with Derive 6.) We observe that our algo- 
rithm produces better results than the ones obtained with the BESSEL.MTH file although the 
proposed algorithm is slower. 


3) Results of the evaluation of the B(z) (See Appendix E). 

The standard BESSEL_Y(k_, z) function cannot give results for the case k_ = 0 and z = 2, but we 
get accurate results with our procedure [5] as it is shown in Appendix E. 

Appendix E 

We have Appendix C loaded and we proceed in order to evaluate the 
Bessel function of the second kind of order 0 up to 10 terms in the summation. 

We take z = 2. (all results are approximated!) 

#20: VECTOR_Y0(10 , 2) 

#21: [0.6366197723, 0.489754084, 0.5119671213, 0.5103024958, 0 . 5103779228 , 

0.5103756229, 0 . 5103756734, 0.5103756726, 0.5103756726, 0.5103756726] 
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We define the Bessel function of the 2nd kind of order 0: 

#22: yO := 0.5103756726 

#23: VECT0R_Y1( 10 , 2) 

#24: [-0.05499896203, -0.1127961327, -0.1066902521, -0.107045282, -0.1070320968, 

-0.1070324379, -0.1070324314, -0.1070324315, -0.1070324315, -0.1070324315] 

We define the Bessel function of the 2nd kind of order 1 : 

#25: yl :i -0. 1070324315 

We evaluate up to 10 terms in the continued fractions expression of H(z=2) with Nmax = 5: 

#26: H C 10 , 5, 2) 

#27: [0, 0.2, 0.2068965517, 0.207070707, 0.2070739549, 0.2070740015, 0.2070740021, 

0.2070740021, 0.2070740021, 0.2070740021, 0.2070740021] 

We define nlim: 

#28: nlim := 6 

We evaluate the Bessel function of the 1st kind of order Nmax: 

#29: 1T0P( 5 , 2) = 0.007039629737 

We evaluate the Bessel function of the 1st kind of order Nmax-1 : 

#30: 1T0PM1(5 , 2) = 0.0339957198 

We evaluate the Bessel functions of the 1st kind of all orders 
below Nmax using the proposed algorithm 

#31: 1NEW( 5 , 2) 

#32: [0.007039629737, 0.0339957198, 0.1289432494, 0.3528340286, 0.5767248078, 

0.2238907791] 

We evaluate the Bessel functions of the 1st kind of all orders 
below Nmax using the BesselFunctions.mth utility file. 

#33: J0LD( 5 , 2) 

#34: [0.007039629755, 0.0339957198, 0.1289432494, 0.3528340286, 0.5767248077, 

0.2238907791] 

7 Conclusions 

We have introduced the Derive version of an algorithm to evaluate the continued fraction representa- 
tion of a rational approximation of a function and an efficient algorithm to evaluate BFs. This code 
illustrates the use of the continued fraction method and gives an accurate way to evaluate the BFs 
although it is slower than the implemented functions in BESSEL.MTH. The proposed method not 
only is very accurate but solves problems as the calculus of, for example, Y 0 (z = 2), that is not possi- 
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ble to obtain the BESSEL_Y(0 , 2) function. Moreover, the aim of this paper is to introduce in an inter- 
active way a complete tested algorithm that can be presented as a tutorial in Bessel functions, claim- 
ing in the stability of the used recurrence relation (i.e. forward recurrence relations for the BFs of 
second kind and backward ones for the BFs of the first kind) and about the properties of convergence 
of the continued fractions approximations. 

Comment in 2008: Times have changed since 2008. The complaint about BESSEL_Y(0,2) is 
obsolete as you can see below. BesselFunctions.mth might be an improved version of 
BESSEL. MTH or it might have been that the authors would have achieved better results 
with the implemented functions by increasing the precision. Josef. 

BESSEI Y(0 , 2) 

0.5103756726 
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Additional comment: 

If you want to know more about Bessel functions and their applications, you can 
find a lot of information in the web. Bessel functions are therefore especially 
important for many problems of wave propagation, static potentials, and so on. 

For example: electromagnetic waves, heat conduction, vibration of membranes, 
diffusion problems, signal processing. Crick and Watson used BFs to develop 
their model of the DNA double helix. 

Three pictures from my Graphic- Workshop held at the DERIVE Days Diisseldorf (DERIVE for DOS). 

Surfaces of Revolution in Parallel- and Central Projection 
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Cubic Splines 

Josef Bohm, Otto Reichel, Leo Klingen 

There are a lot of problems dealing with the task to find a curve connecting more or less given points - 
nodes - in the plane (e.g. road construction, design forms, . . .)• Let us formulate the interpolation prob- 
lem: which continuous and differentiable - not too complex - curve connects n given nodes (x k , y k ) 
with k= 1, . .., n. 


Example: 


k 

Xk 

yk 

Xk+ 1 Xk hk 

yk+ 1 yk jk 

i 

0 

l 

i 

i 

2 

1 

2 

2 

2 

3 

3 

4 

1 

-1.5 

4 

4 

2.5 

1 

-2.5 

5 

5 

0 

1 

1 

6 

6 

1 




The easiest way to find a fitting function would be to find a polynomial of order n- 1 : 

p(x) = a n _ x x n ~ x + a n _ 2 x n ~ 2 + ... + a x x + a 0 . (*) 

The unknown coefficients ao , a\, . . ., a n _\ are the solution of a system of n linear equations: 

[/>(■*! ) = V, , P(x 2 ) = y 2 ,..., P(x n ) = yj. 

Doing this in case of more nodes we might face some disadvantages: 

• a huge system of equation must be solved (could be done using a powerful CAS), 

• it is not very comfortable to work with a polynomial of high order (no problem with a CAS), 

• the curve would show up to n-2 turning points, i.e. the curve could be very oscillating (this 
can not be changed using a CAS). We prefer a smoother curve. 

So let's try another approach tackling the problem. We connect neighbouring points by cubics 

P k (x) = d k + b k {x-x k ) + c k {x-x k ) 2 + d k (x-x k ) 3 

in such a way that they have a node in common (same function values), have the same first derivative 
in this node (then curve is differentiable in all points) and have also the second derivative in common 
(which is responsible that there is no ’’shock” in the curvature). 

(This is a welcome occasion to demonstrate the students the advantage of developing a polynomial on 
a certain position.) 

The easiest way to find the polynomial (*) is applying the FIT-function: 

Let's do this as an intro and then proceed with the cubic splines. 
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We set up 4 linear equations for each of the n - 1 cubics ((1) - (4)). We lack two equations to get a 
unique solution for the system. So we add two meaningful equations. Setting the 2 nd derivative equal 
zero in start and end point is one popular possibility ( natural splines) (equations (5) and (6)). Other 
settings are possible. 

System I: 


Pk( x k)=y k 

k = 1,.. 

.,72-1 

a) 


p k ( x k + i)=yk + i 

k = 1,.. 

,,n — 1 

(2) 


p k '(**+ 1) = p k+ 1 '(**+i) 

k = 1,.. 

. , 72 — 2 

(3) 

!!! 

Pk "(**+i) = A+i "(**+i) 

k = 1,.. 

.,72-2 

(4) 

!!! 

Pi "(*i) = 0 



(5) 


Pn-l "( X « ) = 0 



(6) 



and with 



Pk '(*) = b k + 2c* (x - x t ) + 3</ t (x - 

^) 2 



p k "(x) - 2 c k + 6d k {x-x k ) and h k 



we obtain 



System II: 



(1) 

II 

^3 

k = 1, 

. . . , 72 — 1 

(2) 

a k + b k h k + C k h k 2 + d k h k = JVl 

k = 1, 

. . . , 72 — 1 

(3) 

b k +2 C k h k +3d k h k 2 =b k+l 

k = 1, 

...,72-2 

(4) 

2 C k + bd k b k = 2c k+ \ 

k = 1, 

...,72-2 

(5) 

c\ = 0 



(6) 

2c„_ 1 +6J„_,/?„_ 1 =0 
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For our special problem using the data points System II leads us immediately to System III: 


(1) a x - 1, a 2 - 2, a 3 = 4, a 4 = 2.5, a 5 = 0 

(2) + c x + d x =2 

^ "I - 2Z?2 + 4c?2 + M 2 — 4 
^3 ^3 ^3 “2.5 

a 4 + b 4 + c 4 + d 4 =0 
u 3 -\- b 3 -\- c 3 -\- d 3 =1 

(3) b x + 2c x + 3^ = b 2 
b 2 + 4c 2 + 12J 2 = t> 3 
b 3 + 2c 3 + 3 d 3 — b 4 
b 4 + 2c 4 + 3 d 4 — b 3 

(4) 2c x + 6d x = 2 c 2 
2 c 2 + 12rf 2 - 2 c 3 
2c 3 + 6 d 3 = 2c 4 
2c 4 + 6d 4 = 2 c 5 

(5) c, =0 

(6) 2c 5 + 6d 5 = 0 


We have 4-(6 - 1) = 20 equations. There are 4n - 4 coefficients but only 3n - 4 are really unknown. 
This system can easily be solved - not by hand but by using DERIVE. We obtain the coefficients of 
the p k (x) functions valid for the intervals x k < x < x k +\. Either doing it by hand or using a DERIVE 
function we create the cubics. 



SPLI 1 new.mth 


#1: 

InputMode Word 


#2: 

[nodes := [ ] , coeff :=] 


#3: 

xk VECTOR(nodes , k, DIM(nodes)) 

k, 1 

#4: 

pk := VECTOr[ I coeff 

Li =0 (DIM(xk) - 

•(x - xk ) , k_, DIM(xk) - ll 

l)-i + k_ k_ J 

#5: 

DIM(xk) - 1 

SPL(x) :z I pk . X (xk , 

J=i i j 

x, xk ) 

3 + 1 

#6: 

TAB(incr) := VECTOR([x, SPL(x)], 

x, nodes , nodes , incr) 

1,1 DIM(nodes) , 1 
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#7: Instead of TAB we use now the TABLE-function (see #23 and #26). 
#8: Enter the matrix of the given points and name the matrix nodes. 
#9: coeff is the solution vector - see below. 

#10: pk returns the separated cubic splines and 
#11: SPL(x) gives the complete interpolating function. 


#12: nodes := 


0 1 
1 2 

3 4 

4 2.5 

5 0 

6 1 


lie system: 


#14: [al = 1, a2 = 2, a3 = 4, a4 = 2.5, a5 = 0, cl = 0, al + bl + cl + dl = 2, a2 + 

2 • b2 + 4 • c2 + 8 • d2 =4, a3 + b3 + c3 + d3 = 2.5, a4 + b4 + c4 + d4 = 0, a5 + b5 

+ c5 + d5 = 1, bl + 2 *cl + 3 *dl = b2, b2 + 4*c2 + 12-d2 = b3, b3 + 2*c3 + 3-d3 = 

b4 , b4 + 2 • c4 + 3-d4 = b5, 2-cl + 6-dl = 2-c2, 2*c2 + 12-d2 = 2-c3, 2-c3 + 6-d3 

= 2 «c4, 2 *c4 + 6 • d4 = 2-c5, 2-c5 + 6-d5 = 0] 
coeff := (S01UTI0NS( [al =1, a2 = 2, a3 = 4, a4 = 2.5, a5 = 0, cl = 0, al + bl + cl 

#15: 


+ dl = 2, a2 + 2 «b2 + 4-c2 + 8-d2 = 4, a3 + b3 + c3 + d3 = 2.5, a4 + b4 + c4 + 


d4 =0, a5 + b5 + c5 + d5 = 1, bl + 2-cl + 3-dl = b2, b2 + 4-c2 + 12-d2 = b3, b3 


+ 2 • c3 + 3 *d3 = b4, b4 + 2*c4 + 3-d4 = b5, 2*cl + 6-dl = 2-c2, 2-c2 + 12-d2 = 


2 *c3 , 2 *c3 + 6*d3 = 2-c4, 2-c4 + 6-d4 = 2-c5, 2-c5 + 6-d5 = 0], [al, a2 , a3 , a4, 
a5 , bl, b2, b3, b4, b5, cl, c2, c3, c4, c5, dl, d2, d3, d4, d5])) 


#16: coeff := 


1, 2, 4, , 0, 

2 


45 111 


38 38 ' 38 


19 


33 

24 

6 

38 

19 

19 

0, 

26 

37 

19 

38 



51 


19 


18 


19 


0, 


1 

15 

38 


45 


38 


#17: pk 


# 18 : 


3 3 2 

5-x 33- x 10-x - 45-x + 12-x - 53 

+ + 1 , - 


38 


38 


38 


45-x - 258-x + 217 


38 


3 2 2 

52-x - 669 • x + 2754* x - 3545 (5 - x).(37-x - 481-x + 1516) 


38 


38 
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#19: SPL(x) 


2 2 
(x - 5). (37-x - 481- x + 1516) • SIGN(x - 6) 89-(x - 10-x + 25). |x - 5| 


2 2 
13* (x - 8-x + 16). |x - 4 | 5 • (x - 6-x + 9) • | x - 3| 


2 3 

15- (x - 2-x + 1) * | x - l| (5 • x + 33-x + 38).SIGN(x) 


#21: TAB(0 . 5) 


0 +0.5 + 0.5 


0.5 1.450657894 


#23: TABLE(SPL(x) , x, 0, 6 P 0.5) 


0 +0.5 + 0.5 


0.5 1.450657894 


1.5 2.697368421 


1.5 2.697368421 


2 3.394736842 


2.5 3.894736842 


#22: 3 


3.5 3.546052631 


2 3.394736842 


2.5 3.894736842 


#24: 3 


3.5 3.546052631 


4.5 1.032894736 


4.5 1.032894736 


5.5 0.1348684210 


6 +0.5 + 0.5 


5.5 0.1348684210 


6 +0.5 + 0.5 


The function values in the nodes are strange - this is a consequence of using the CHI-function. We 
could overcome this problem by using an IF-construct to define the piecewise defined function. 

It is interesting that we can close the gaps by redefining the function: 


#25: SPL_(x) :z 


(x - 5) • (37 • x - 481- x + 1516) * SIGN(x - 6) 


2 2 
89- (x - 10- x + 25) - | x - 5 | 13-(x - 8-x + 16). |x - 4| 


2 2 
5-(x - 6-x + 9) - | x - 3 | 15. (x - 2-x + 1) • | x - l| 


(5-x + 33-x + 38) « SIGN(x) 


76 
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#26: TABLE(SPL_(x) , x p 0 P 6, 0.5) 

Mi 0 ± 0 . 5 + 0 . 5 


#27: 


0.5 1.450657894 

1 2 

1.5 2.697368421 

2 3.394736842 

2.5 3.894736842 

3 4 

3.5 3.546052631 

4 2.5 

4.5 1.032894736 

5 0 

5.5 0.1348684210 

6 ±0.5 + 0.5 


#28: The polynomial of order 5 through the nodes: 


#29: 


#30: 


15 4 3 2 I 

FIT(Lx, a-x + b-x + c-x + d-x + e-x + fj p nodes) 

5 4 3 2 

x 5 • x x 29- x 17- x 

- + + + + 1 

80 48 16 48 40 


#31: 

#32: 


You can load this file as a utility file! 


rABLE(SPL_(x) , x, 0, 6, 0.01)1 


The last expressions is the base for plotting the spline as a thick line (Points connected and size large). 


But you might prefer another way. The next procedure is creating the system of equations, too: 

SPLI_2_new . mth and SPLI_2_new . dfw 

#1: d(a) := DIM(a) - 1 

[ x_(a, k) := a , y_(a, k) := a 1 
#2: L k , 1 k , 2 J 

#3: h(a, k) := x_(a, k + 1) - x_(a, k) 

#4: m0(a) :r VECTOR(VECTOR(IF(k = 4-j - 3, 1, 0), k, 4-d(a)), j, d(a)) 

#5: ml(a) :r VECTOR(VECTOR(IF(k = 4-j - 3, 1, IF(k = 4-j - 2, h(a, j), IF(k = 4-j - 1, 

2 3 

h(a, j) , IF(k r 4-j, h(a, j) , 0)))), k, 4-d(a)), j, d(a)) 

#6: m2(a) :r VECTOR(VECTOR(IF(k = 4-j - 2, 1, IF(k r 4-j - 1, 2-h(a, j), IF(k r 4-j, 

2 

3-h(a, j) , IF(k = 4-j + 2, -1, 0)))), k, 4-d(a)), j, d(a) - 1) 
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8 7 6 5 4 3 2 

#23: FITCLx, a-x + b-x + c_-x + d_-x +e-x + f-x +g-x + h_-x + ij, pts) 

-5 8 7 6 5 4 

#24: 6.109773165-10 -x - 0. 0008684255363 -x - 0. 001684734074 -x + 0. 05405121871- x - 0. 03029723448- x - 

3 2 

0. 9550514929- x + 0. 5513712830 -x + 5 . 510965706- x + 0.3714525808 



You might guess which one of the two curves is the spline and which one is the polynomial? 

a represents the matrix of the points, system (a) gives the left hand side matrix of System III from 
above, c0(a) gives the right hand side of System III (the constants), coef f (a) returns the solutions, 
spl (a) shows the cubics in a vector and finally spl ines(a) creates the composed function valid for 
the whole interval. 


With a bit of manipulation skills System III can be transformed into a special form. This could be 
done with the students and presents a fine example that even in times of a CAS some manipulation 
skills are very useful. 

We would like to eliminate all the variables in System II except the c' s. In the first two steps we 
eliminate a k and d k using (1) and (4) (page 28), then we set y k+x -y k =j k . Now the two equations (2) 
and (3) are remaining (containing b k and d k ): 


(2) b k h k + c k h k 2 + C 




3 


■ j k l solve for b k b k =... (*) 


K=- 


(3) b k + 2 c k h k + (c t+1 - c k ) h k = b k+] ; solve for b k 

We compare the two expressions for the b k and solve the resulting equation for b k+ \\ 

b =^ + c h - ° k+] ~° k h 

u k + i 7 ^ c k + i n k ~ n k • 

K 3 

Decreasing the index leads to a new expression for b k . 


b = dkzL + c h __ 

u k 7 ^ c k n k - 1 

K-x 


If we now equate this expression with the expression (*) for b k we obtain System IV: 


c k- 1 K-\ + ^ c k {h k + h k _ x ) + c k+l h k — 3 

= C n =0 


f ^ 

J k J k- 1 

h h 

\ n k n k - 1 


\ k = 2,...,n-l 
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= 0 

(k = 2) 

= -7.5 

II 

UJ 

= -3 

(k = 4) 

= 10.5 

5? 

ii 


And in our case System V at last is hiding the spline problem's solution: 

6c 2 + 2 c 3 
2 c 2 + 6c 3 + c 4 

c 3 + 4 c 4 + c 5 
c 4 + 4c 5 

For solving this symmetric "tridiagonal" system exists a fast algorithm even for large systems and 
having found the cr coefficients it is not difficult to find the values for the b k , the a k and d k . (See ref- 
erences). 

For the real DERIVIANs it might be nice to perform the manipulations from above using the CAS. Stu- 
dents could check their competence in manually manipulations and DERIVE manipulations as well: 

I substitute b k = b_: 


#1: SOLVE 


c - c 

2 k+1 k 2 

b_-h + c -h + h = ] , b. 

k k k 3 k k J 

2 2 
2-c-h + c -h - 3 - j 

k k k+lk k 


#2: b_ = - 


3-h 


SOLVE (b_ + 2-c -h + (c -c)-h = b , b_) 

#3: kk k+1 kk k+1 

b_ = b - c -h - c -h 

#4: k+1 k+lk kk 

2 2 

2-c-h +c -h - 3 - j 

k k k+lk k 

#5: - = b -c -h-c-h 

3-h k+1 k+lk kk 

k 

I substitute b k+1 = b . (In Simplify Menu: Subexpression Substitution) 

2 2 

2-c -h + c -h - 3-j 

k k k+lk k 

#6: - = b - c -h - c -h 

3-h k + 1 k k k 


#7: SOLVE 


2 2 
2-c -h + c -h - 3-] 
k k k+lk k 

3-h 


b_ - c -h - c -h , b_ 
k+lk k k 


#3: b = 


2 2 

2-c -h +c-h +3-] 
k+lk k k k 

3-h 


I substitute b -> b_ and k -» k - 1 and then simplify expression #10: 


#9: b_ 


2 2 

2-c -h +c -h +3-] 

(k - 1) + 1 k - 1 k-lk-1 k - 1 


3-h 


k - 1 
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#10: b_ 


2 2 

c ■ h + 2 - c - h + 3 - ] 

k-lk-1 kk-1 k - 1 

3-h 


k - 1 

2 2 2 2 
2-c-h + c -h - 3 - j c ■ h + 2 - c - h + 3 - ] 

kk k+lk k k-lk-1 kk-1 k-1 


# 11 : - 


3-h 


3-h 


k-1 


2 2 2 2 
2-c-h + c -h -3-] c -h + 2-c -h + 3 - j 

kk k+lk k k-lk-1 kk-1 k-1 


#12: sy5lV := - 


3-h 


3-h 


k k-1 

Expression #13 creates the system, h and j are vectors obtained from the table page 27. 

#13: VECTOR(sysIV, k, 2, 5) 

#14: [It := [1, 2, 1, 1, 1], j :: [1, 2, -1.5, -2.5, 1]] 


#15: 


2-c+c 8-C+4-C+15 

2 1 3 2 

c + 2-c = - , c + 2-c = - , c + 2-c = - 2-c 

32 243 2 544 


4-c + 2-c - 21 

5 4 


c - 3, c + 2-c 
3 6 5 

I substitute Ci = c 6 = 0: 


#17: 


#18: 


2-c + 0 8-C+4-C+15 

2 3 2 

c + 2-c = , c + 2-c = - , c + 2-c = - 2-c - 

32 243 2 544 

4-c + 2-c - 21 

5 4 

c - 3, 0 + 2-c = 

3 5 2 

8-c + 4-c + 15 

3 2 

c + 2-c = - c , c + 2-c = , c + 2-c = - 2-c - c - 3, 

32243 2 54 43 

4-c + 2-c - 21 

5 4 


5 4 

Unfortunately we can not recognize the nice "tridiagonal" form in expression #18. But we try to solve 
this system for the Cj. I replace the Cj by cj and solve the system: 

#19: InputMode := Word 


#20 : SOLVE 


8 • c_3 + 4 • c_2 + 15 

c_3 + 2 -c_2 = -c_2 , c_4 + 2-c_3 = , c_5 + 2-c_4 = 


4-c_5 + 2 -c_4 - 21 


2 -c_4 - c_3 - 3, c_5 = 


# 21 : 


15 45 45 

c_2 = a c_3 = - a c_4 = a c_5 = 

38 38 38 38 


, [c_2 , c_3 , c_4, c_5] 
111 


Comparing this with the solution of SYSTEM V or with the respective values in the coeff-vector in ex- 
pression #16 page 30 we can be proud that we did a good job manually and "CASially" as well. 


I continue with transferring the said algorithm into a DERIVE code. 
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# 1 : 

#2: 

#3: 

#4: 


file spl i_3_new. mth 


nodes := 
h(k) := nodes 


- nodes 
k + 1,1 k, 1 


f(k) := 2 • (h(k + 1) + h(k)) 


#5: m(k) := 


nodes - nodes 

k + 1,2 k , 2 


# 6 : 

#7: 

#8: 

#9: 

#10: 


h(k) 

g(k) := 3 • (m(k + 1) - m(k)) 

a(n) := 

If n = 1 

f(l) 

f(n) - h(n) A 2/a(n - 1) 

aO :: APPEND(VECTOR (nodes , i, DIM(nodes) - 1), [nodes 

i .2 L 

a_val := VECTOR(a(i), i, DIM(nodes) - 2) 

b(n) := 

If n = 1 

g(D 

g(n) - b(n - l)-h(n)/a(n - 1) 


1) 

DIM(nodes),2j 


#11: b_val := VECTOR(b(i ) , i, DIM(nodes) - 2) 

c_aux(n) := 

If n = 1 

#12: b(DIM(nodes) - 2)/a(DIM(nodes) - 2) 

(b(DIM(nodes) - n - 1) - h(DIM(nodes) - n)-c_aux(n - l))/a(DIM(nodes) - n - 1) 

#13: a2 := APPEND( [0 ] , VECTOR (c_aux(j) , j, DIM(nodes) - 2, 1, -1), [0]) 

r a2 - a2 

] + 1 j 


#14: a3 := VECTOR 


#15: al := VECTOR 


j, DIM( nodes) - 1 


3-h(j) 

r aO - aO h(])-(a2 + 2*a2 ) 

j + 1 j j + 1 j 


, ], DIM( nodes) - 1 


#16: splines := VECTOR 


h(]) 


aO + al *(x - nodes ) + a2 *(x - nodes ) + a3 *(x 

V j j j,l j j.l J 


nodes ) , j, DIM(nodes) - 1 

]M 

DIM(nodes) - 1 

#17: spl(x) := I splines *x(nodes , x, nodes ) 

i=l i i , 1 i + 1, 1 

0 1 

1 2 

3 4 

4 2.5 

5 0 


#18: nodes := 


6 1 
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#19: splines 


3 2 

33 -x 10- x - 45 -x + 12- x - 53 

+ 1 , 

38 38 


2 

45 -X - 258 ■ > 

c + 217 

3 2 

52 -X - 669 -X + 2754 ■ X - 3545 

38 

1 

38 


(5 - x) • (37 -x - 481- x + 1516) 
38 


#20: spl(x) = 


(x - 5) • (37 *x - 481-x + 1516) • SIGN(x - 6) 


76 


3 2 

89- (x - 15-x + 75-x - 125) • SIGN(x - 5) 


76 


3 2 

13- (x - 12-x + 48-x - 64) • SIGN(x - 4) 


19 


3 2 3 2 

5 • (x - 9 • x + 27- x - 27) «SIGN(x - 3) 15-(x - 3-x + 3-x - l).SIGN(x - 1) 


38 


76 


(5-x + 33-x + 38) • SIGN(x) 


76 


It is very informative to plot the single cubics in different colours - together with the resulting spline. 
One can see the very smooth transition from one cubic to the next one. 



Starting with DERIVE 5 we could program and so we are able to collect the whole procedure 
into one program. As an extra you can find Don Phillip's program among the accompanying 
files. Another program can be found in my book "Programmieren in Derive" (bk-teachware). 
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At last Yd like to add that I like the "cubic splines”, because the "Reverse Curve Discussion" (in Aus- 
tria: Umgekehrte Kurvendiskussion, in Germany "Steckbriefaufgaben" = "wanted circular problems") 
are a flxpoint in our school mathematics and I for my person have some problems answering questions 
about the importance of these problems in real life. I can assure you that students also are enthusiastic 
with the splines. 

In DNL#19 I'll show how Otto Reichel from St. Polten and Leo Klingen from Bonn are dealing with 
this issue. As a special gift you will find Mr Klingen's function to produce paramtetric splines. (Do 
you remember the Easter Bunny from last year?) 

At the DERIVE Days Diisseldorf Gunter Scheu promised to send a paper on Bezier Curves. And he 
really did. Thank you, Gunter, you are very reliable. So weTl see how this interpolation technique is 
working soon. 

I found a nice website with an applet for calculating and presentating cubic splines: 
http://www.amdt-bruenner.de/mathe/scripts/kubspline.htm 



Berechnen 


Krummung in Anfangs- und Endpunkt festlegen : 
Erster Punkt: S"(xo)= 0 


3.5 


Letzter Punkt: S"(x n )= 


Interpolieren 


Koeffizienten: 


B ruche (exakt) 


S(x) = 


3 ,546052631579 


x au3 [0; 1] a 

50 (x) = 5/ 3 3 ■ x A 3 + 33/38-x + 1 

= 5/33 ■ x A 3 + 33/38-x + 1 
x aus [ 1 ; 3 ] 

51 (x) = -5/ 19 ■ (x-1) A 3 + 15/38- (x-1) A 2 + 24/19- (x-1) + 2 

= -5/ 19 ■ x A 3 + 45/3 3 ■ x A 2 - 6/19-x + 53/38 


Nullstellen: { 5; 5,369943794749 } 

Extremstellen: lok. Minimum: { 5,178005063473213 } - lok. Maximum: < 2,860147050873544 } 
Wendestellen: { 0; 1,5; 4,288461538461538; 6 } 


Flachenintegral 

5,36994379 

J S(x) dx = -0,020089445825 


Volumenintegral 

5,36994379 

71 J S(X)2 dx = 0,004116102745 


Kurvenintegral (approximiert) 

6 

J V(l+S'(x)2) dx 10,529 
o 

r^i automatisch r-z- 

^ berechnen Genauigkeit: 10 LA 
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You may remember my "DERIVE program" for investigating a curve from DNL#15. Otto Rei- 
chel has written a related DERIVE tool for the "reverse" discussion of a curve. This is the 
exact translation word by word how we - teachers and students in Austria - often name this 
kind of problems. In Germany they have another name "Wanted - Problems". 

Let 's take a typical problem: The graph of a polynomial function of order 4 passes two points 
P(2|0) and Q(-1|0). The origin is an inflection point with 1 as slope of the tangent. Find the 
equation of the function. 

In my opinion Otto's file does not need any further explanation. The first method which he 
uses is very common and not new, but the second one provides a skilled function. DERIVE 
"looks" for the order of the polynomial automatically and within one step you will obtain the 
result. This is really a fine tool for teachers and a nice programming technique. Josef 


The "Reverse" Discussion of a Curve 

Otto Reichel, St. Polten, Austria 


z 

#1: fZ_(x) := a x + b-x + c 

#Z : fZ_l(x) := Z ■ a ■ x + b 

3 Z 

#3: f3_G(x) := ax + b x + c-x + d 

Z 

#4: f3_i(x) :=3 a x + Zbx + c + d 

#5: Read as: Polynomial of degree 3 - 1st derivative 

#6: f3_Z(x) := 6- a x + Z-b 

4 3 Z 

#7: f4_0(x) := ax + b-x + ox + d x + e 

3 Z 

#3: f4_l(x) := 4-a-x + 3-b-x + Z-ox + d 

Z 

#9: f4_Z(x) :: IZ-a-x + G-b-x + Z-c 

#10: Example 1: 

#11: The Graph has two zeros at xl=G and xZ=Z , P(l|-1) is an inflection point 
#1Z : with the slope of the tangent = -Z . 

#13: SOLVE ( [ f 4_0 (0 ) = 0, f4_0(Z) = 0, f4_0(l) = -1, f4_Z(l) = 0, f4_l(l> = -Z ] , [a, b, c a d a e]) 

#14: [a - 1 a b = -Z ac =0 a d = 0 Ae =0] 

4 3 

#15: x - Z-x 



#19: (SOLUTIONS([f4_O(0) = 0, f4_0(2) = 0, f4_0(l) = -1, f4_0"(l) = 0, 


f4_0'(l) = -2], [a, 


b, c, d. 



3 

X f 


2 

X f X | 


J 


4 3 

#20: x * 2'X 


#21: The general approach is following. 


#22: 

h(n, x_, k) := VECTOR^ 

trzF - "■ - 1 ) 

#23: 

' h(4, x_, 0) ■ 
h(4, x_, 1) 


[. 

r 4 3 2 i 

Lx_ , x_ , x_ , X_, 1J 

3 2 1 
4*x_ , 3*x_ , 2*x_, 1, Oj 


. h(4, x_, 2) . 



[ 12 • x_ , 6-x_, 2, 0, o] 


mat(in) := VECT0R( lim h(DIM(in) - 1, x_, in ), i, DIM(in)) 
#24: x_-»inii|2 i,l 


-1 

#25: coeff(in) := mat(in) *in COL 3 

#26: func(in) := coeff(in)»h(DIM(in) - 1, x, 0) 


#27: Example 1: 


#28 : matexl 


0 0 0 

0 2 0 

2 10 

0 1-1 


L 1 1 -2 J 


4 3 

#29: func(matexl) = x - 2-x 

#30 : Example 2 : 

#31: The graph of a polynomial of degree 4 (quart ic) which is symmetric wrt y-axis 

#32: has an inflection point in I(-2|3) with slope 4. 


#33 : matex2 := 


0- 2 3 
2-2 0 

1- 2 4 

0 2 3 


2 2 0 
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4 2 

x 3 ■ x 

#34: func(matexZ) = + 8 

16 Z 


#35: func 


0 -Z 3 

Z -Z 0 

1 -Z 4 

0 Z 3 


4 Z 

x 3 ■ x 

- + 8 

16 Z 


Z Z 0 

1 Z -4 . 


#36 : Example 3 : 

#37: Find the polynomial function with minimum order which shows 

#33: three turning points in T1(0 |Z), TZ(1 | —3 ) and T3(4|4>. 


0 0 z 
10 0 


#39: func 


0 1-3 

110 


0 4 4 


14 0 



On the next page you can find a ready made tool for finding the polynomial functions for the 
Tl-devices and NspireCAS. The Voyage 200/TI-92 function is the same as the code pre- 
sented in the Tl-Nspire function which follows. There is only one difference: instead of a -» b 
for storing a as b on the TI-92/V 200 one has to enter b := a with Tl-Nspire. Josef 
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r-ff 

fet y f::t y fht y fe y fs^ y 
H 1 gebra |Calc lot-hen |PrgnI0|c lean Up| | 

■ fnc 

■0 -2 3 ‘ 

2-2 0 

1 -2 4 

0 2 3 

2 2 0 

.1 2 -4. 

x 4 3x2 +3 

IS 2 + 8 

^ ..3 

...2.4][G.2.3][2.2.0]E1.2.-4]]> 

MAIN RAD EKACT FUNC 30/30 


f " 1 gebralc"lcIotherIprgm I ole 1 e"n Up ! 1 

u 

ui 

.1 1 - 2 .. 

'0 0 0 " 

0 2 0 

1 

■ fnc 

2 1 0 
0 1 "I 

.1 1 - 2 .. 

x 4 - 2 x 3 

... 0,2 

:,Q;2,i, 

,0;G-l--l;l-l.-2]> 

MAIN 

RAD EKACT FUNC 30/30 



fnc 


0 2 
0 0 
1 -3 
1 0 
4 4 
4 0 


313*x 

432 


64-x 4 3175-x 3 

+ 

9 144 


4463'x 

216 


-+2 


D Graph passes (— 2|3); turning point in inflection point in (3,- 

slope at x=-4: — 

2 


i); 


fnc 


1 -4 - 


44767'X 1961 lvx 3G791-x J 181221-x"' 136949-x 3144 


458880 


152960 


11472 


-+- 


38240 


-+- 


14340 239 


3/99 


0/99 


fnc 


10/10 


Define LibPub fnc( m_)= 

Func 

Local mfnjjtXp 

n : =dim(m_j[l] : m 1 : =newMat(n, n+ 1 ) 

For i.l.n 
For Jjljfl 

— r-pr(x"“-')|x=m_[i J 2] 

cbT- L'UJ 

EndFor 
EndFor 

xp \= 1 ist. ► mat(seq (x* k, n - 1 , G, - 1 )) 
dotp((rref(ff? j)) T [tf+l] J x j p) 

EndFunc 
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Titbits from Algebra and Number Theory 

by Johann Wiesenbauer, Vienna 

In the last article of this series I reported about some minor bugs in version 3.04 of DERIVE. In the 
meantime all of them have been removed, as I am happy to tell you. Therefore lucky owners of 
DERIVE 3.05f (or later) should use the updated versions of the following routines in that article 
(thereby also correcting some misprints of mine): 

p(n) := Ff(IF(PRIME(k_) , -1, 0), k_, FACTORS (FACTOR(n))) 

SQUAREFREE(n) := IF(n(n) A 2 = 1, true, false) 

5(n) := n(ITERATE(IF(PRIME(f_ A (l/k_)) , k_, k_, k_ +1, k_ + 1) , 
k_, 1) + 1, f_, FACTORS (FACTOR(n))) 

cpCn) := n-nCl - l/f_ A (l/ITERATE(IF(PRIME(f_ A (l/k_)) , k_, k_ + 1, k_ + 1) , 
k_, FACTORS (FACTOR(n))) 

LEGENDRE(a, p) := IF(- PRIME(p) v p = 2, "Input error!", 

M0DS(a A ((p - l)/2) , p)) 

PEPINTEST(m) := IF(m = 0, true, IF(MODS(3 A 2 A (2 A m - 1), 2 A 2 A m + 1) = -1, 
true, false)) 

In particular, as you may conclude from the last two examples, the MODS-function uses the power- 
mod-algorithm as well now, when it comes to calculating powers. Great! 

Referring to the contribution of Bernhard Wadsack in DNL#17, 1 would like to point out that DERIVE 
cannot only display huge Mersenne primes but also prove that they are prime! The basic tool to do this 
is the so-called Lucas-Lehmer test which can be stated as follows: 

If p is an odd prime, then the Mersenne number M p = 2 P - 1 is prime if and only 
if s p i = 0 mod M p , where the sequence s i, s 2 , s 2 , ... is defined by 

s \ =4, s k = sl_ x - 2 modM p . 

Its implementation in DERIVE could look like this: 

LLTEST(p) := IF(p = 2 v PRIME(p) , "p is supposed to be an odd prime", 

IF (ITERATE (MOD (s_ A 2 - 2 , 2 A p - 1) , s_, 4 , p - 2) = 0 , true, false)) 

As a first test we will check an assertion of Mersenne in 1644, who claimed that the numbers 

M p = 2 P - 1 are prime for numbers p < 257 if and only if p is one of the following primes: 2, 3, 5, 7, 

13, 17, 19,31,67, 127, 257. 

For the sake of convenience, we write a little routine that checks all Mersenne numbers 
M p = 2 P - 1 with p < s for primality: 

MERSENNE_LIST(s) := SELECT(LLTEST(p_) , p_, SELECT(PRIME(k_) , k_, s)) 
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The following line was executed in 8.1 seconds (in 1995, 0.172 sec in 2008): 
MERSENNE_LIST(257) = [3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127] 


As it turned out, Mersenne made 5 errors by including p = 67, 257 and leaving out 61, 89, 107. 


Here are some more examples together with the respective calculation times: 


LLTESTC521) = true (1.0 sec) 
LLTESTC607) = true (1.5 sec) 


LLTESTC1279) 

= true 

(6.7 sec) 

LLTESTC2203) 

= true 

(28.5 sec) 

LLTESTC2281) 

= true 

(31.7 sec) 

LLTESTC3217) 

= true 

(83.1 sec) 

LLTESTC4253) 

= true 

(183.4 sec) 

LLTESTC4423) 

= true 

(207.8 sec) 

LLTESTC9689) 

= true 

(2092.6 sec) 


(0.031 with DERIVE 6.10) 
(0.031 with DERIVE 6.10) 
(0.11 with DERIVE 6.10) 
(0.25 with DERIVE 6.10) 
(0.27 with DERIVE 6.10) 
(0.50 with DERIVE 6.10) 
(0.97 with DERIVE 6.10) 
(1.06 with DERIVE 6.10) 
(9.94 with DERIVE 6.10) 


Though the calculation times are increasing very rapidly, the primality testing of Mersenne numbers 
with several thousand digits is well within reach! The reader might try some more exponents of the 
following list, yielding Mersenne primes, namely 


p = 9941, 1 1213, 19937, 21701, 23209, 44497, 86243, 1 10503 
132049, 216091, 756839, 859433. 


(Well, to be honest, I wouldn't recommend taking numbers of the second half of this list!) 

We turn to a more interesting question, namely: How can we find actual factors of a Mersenne num- 
ber which has proven to be composed by the Lucas-Lehmer test? In general, this is a very difficult 
question. 

Let's take as an example the case p = 67 from above where Mersenne was mistaken in believing that 
it yields a Mersenne prime. Applying factor to M 6 7 - 1 yields in 1.8 sec (!) the following decomposi- 
tion 


FACT0R(2 a 67 - 1) = 193707721-761838257287 (now in 0.078 sec!) 

(By the way, F. N. Cole gave a talk about this factorization at a meeting of the AMS in 1903; accord- 
ing to his own words it took him the "Sundays of three years" to find it.) The following calculations 
show that both factors are = 1 mod 2 p and = +1 mod 8, which is true in general: 

3 3 

FACTOR (193 707721 - 1) = 2 -3 -5- 67- 2677 

2 

FACTOR (76183 32 57237 - 1) =2-3 -29- 67- 2551- 3539 


M0DS(193707721, 3) a 1 
M0DSC76133S257287, 8) = -1 

This fact could be used to speed up trial division. In case of the first prime factor q we also notice that 
q - 1 splits up into small prime factors apart from a single bigger one. There is a factorization method 
introduced by Pollard called (p - l)-method which takes advantage exactly of this fact. 
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Suppose that for a given number n and some boundaries si, s2 all prime powers contained in n are 
less than s 1 apart from a single prime which must be less than s2. By forming appropriate powers a 
of an integer a > 1 and calculating the GCD(a r - 1 ,ri) eventually a factor of n will be found. The in- 
terested reader will have no difficulty in finding out the details by looking into the following 
DERIVE-routine. 

PMINUSl(n, a, si, s2) := (ITERATE(IF(c_ > s2 OR e_ > 1, [a_, b_, c_, d_, e_] , 
[M0D(a_-M0D(b_ A (d_ - c_) , n) , n) , b_, d_, NEXT_PRIME(d_) , GCD(a_ - 1, n)]), 
[a_, b_, c_ , d_, e_] , ITERATE(IF(c_ > si OR e_ > 1, [a_, b_, c_, d_, e_] , 

[IF Cd A 2 > si, MOD(a_ A d_, n) , M0D(a_ A d_ A FL00R(LN(sl) , LN(d_)) , n)), a_, d_, 

NEXT_PRIME(d_) , GCD(a_ - 1, n)]), [a_, b_, c_, d_, e_] , 

[a, 1, 1, 2, GCD(a, n)])))|5 

If n is of the form ( 2'" - 1 )/f where /is some already known factor, we should start with a p instead of 
a. This is an adoption of the (p - 1 /method taking this into account: 

MPMINUSlCp, f, a, si, s2) := PMINUSl((2 A p - l)/f, a A p, si, s2) 

Now let's try if it works: 

MPMINUS1C67, 1, 30, 50, 10000) = 193707721 (2.7 sec) (O.llsecnow) 

Terrific! This execution time comes even close to the time of the built-in factor-routine which is al- 
ready amazingly good! 

The (p - l)-method works also for Fermat numbers F m = 2 +1 which have turned out to be com- 

posed by Pepin's test. The other parameters in the following routine FPMINUS1 have the same mean- 
ing as in MPMINUS1: 


r 


FPMINUSl(m , f, a, si, s2) := PMINUS1 




m 

2 

2 + 1 


f 


m + 2 

2 

a , si, s2 


FPMINUS1C10, 1, 3, 40, 100) = 6487031809 (0.031 sec) 

14 2 

FACTOR (6487031809 - 1) = 2 -3 -29 -37- 41 


For Fermat numbers F m = 2 2 +1 every factor is of the form k 2 m+2 + 1 for some k. The following 
routines calculate factors of F m and the corresponding k. 

LDF(m, si, s2) := ITERATE(IF(MOD(2 A 2 A m , t_) = t_ - 1, t_, IF(t_ > s2, 1, 
t_ + 2 A (m + 2))), t_, IF(M0D(sl, 2 A (m + 2)) = 1, si, 

(FL00R(sl - 1, 2 A (m + 2)) + l).2 A (m + 2) + 1)) 


KLDF(m, si, s2) := (LDF(m, sl«2 A (m + 2) + 1, s2*2 A (m + 2) + 1) - l)/2 A (m + 2) 
KLDFC1945, 1, 10) = 5 (0.48 sec) 


The Fermat number F m5 can never be seen in full length because the number of particles in our uni- 
verse would not suffice to print it, as Coxeter pointed out. Even so, we have just calculated the factor 
5 • 2 1947 + 1 of it! More about Fermat numbers and their factors another time! 


